Data
# ********************************** Data **************************************
hcog <- read.csv("~/Psychometrics Conference/2024/Simulation WG/Data/HRS_Merge3.csv")
names(hcog) <- tolower(names(hcog))
# MMSE
hcog <- hcog %>% mutate(
ornttime_mmse = r1orient_time,
orntplace_mmse = r1orient_place,
imrec_mmse = r1imrc3,
delrec_mmse = r1dlrc3,
spellbck_mmse = backspel,
naming_mmse = r1object,
phrase_mmse = case_when(
!is.na(r1repeat) ~ case_when(
r1repeat == 1 ~ 1,
r1repeat == 5 ~ 0
)
),
commfoll_mmse = case_when(
!is.na(r1combfol) ~ case_when(
r1combfol %in% c(1,2) ~ 1,
r1combfol == 5 ~ 0
)
),
writesent_mmse = case_when(
!is.na(r1senten) ~ case_when(
r1senten %in% c(1,2) ~ 1,
r1senten == 5 ~ 0
)
),
draw_mmse = case_when(
!is.na(r1draw) ~ case_when(
r1draw == 1 ~ 1,
r1draw == 5 ~ 0
)
),
comm3step_mmse = follow3_step,
mmsetot_mmse = r1mmse_score
)
# HRS TICS
hcog <- hcog %>% mutate(
wlimrc_tics = pd174,
wldlrc_tics = pd184,
# count_tics = pd170,
animfl_tics = pd196,
numser_tics = pd216,
serial7_tics = ser7sum,
ornttime_tics = case_when(
!is.na(pd151) ~ case_when(
pd151 == 1 ~ 1,
pd151 == 5 ~ 0
)
) +
case_when(
!is.na(pd152) ~ case_when(
pd152 == 1 ~ 1,
pd152 == 5 ~ 0
)
) +
case_when(
!is.na(pd153) ~ case_when(
pd153 == 1 ~ 1,
pd153 == 5 ~ 0
)
) +
case_when(
!is.na(pd154) ~ case_when(
pd154 == 1 ~ 1,
pd154 == 5 ~ 0
)
),
naming_tics =
case_when(
!is.na(pd155) ~ case_when(
pd155 == 1 ~ 1,
pd155 == 5 ~ 0
)
) +
case_when(
!is.na(pd156) ~ case_when(
pd156 == 1 ~ 1,
pd156 == 5 ~ 0
)
) +
case_when(
!is.na(pd157) ~ case_when(
pd157 == 1 ~ 1,
pd157 == 5 ~ 0
)
) +
case_when(
!is.na(pd158) ~ case_when(
pd158 == 1 ~ 1,
pd158 == 5 ~ 0
)
),
cntbck_tics = pd170 - ornttime_tics - naming_tics
)
# HCAP
hcog <- hcog %>% mutate(
naming_hcap =
case_when(
!is.na(r1scis) ~ case_when(
r1scis == 1 ~ 1,
r1scis == 5 ~ 0
)
) +
case_when(
!is.na(r1cactus) ~ case_when(
r1cactus == 1 ~ 1,
r1cactus == 5 ~ 0
)
) +
case_when(
!is.na(r1pres) ~ case_when(
r1pres == 1 ~ 1,
r1pres == 5 ~ 0
)
) +
case_when(
!is.na(r1elbow) ~ case_when(
r1elbow == 1 ~ 1,
r1elbow == 5 ~ 0
)
) +
case_when(
!is.na(r1hammer) ~ case_when(
r1hammer == 1 ~ 1,
r1hammer == 5 ~ 0
)
),
ceradtot_hcap = r1word_total,
ceraddel_hcap = r1word_dscore,
ceradrcg_hcap = r1wlrec_totscore,
animfl_hcap = r1verbal_score,
bmim_hcap = r1bm_immscore,
bmdel_hcap = r1bm_delscore,
lmim_hcap = r1lmb_immscore,
lmdel_hcap = r1lmb_delscore,
lmrcg_hcap = r1lmb_recoscore,
conprxim_hcap = r1cp_score,
conprxdel_hcap = r1cpdel_score,
dsmt_hcap = r1dig_score,
numser_hcap = r1ns_score,
ravens_hcap = r1rv_score,
trailsa_hcap = r1tma_score,
trailsb_hcap = r1tmb_score,
spatial_hcap =
case_when(
!is.na(r1store) ~ case_when(
r1store == 1 ~ 1,
r1store == 5 ~ 0
)
) +
case_when(
!is.na(r1point) ~ case_when(
r1point == 1 ~ 1,
r1point == 5 ~ 0
)
)
)
mmse_nms <- c("ornttime_mmse","orntplace_mmse","imrec_mmse","delrec_mmse",
"spellbck_mmse","naming_mmse","phrase_mmse","commfoll_mmse","writesent_mmse",
"draw_mmse","comm3step_mmse")
tics_nms <- c("wlimrc_tics","wldlrc_tics","animfl_tics",
"numser_tics","serial7_tics","ornttime_tics","naming_tics","cntbck_tics")
hcap_nms <- c("naming_hcap","ceradtot_hcap","ceraddel_hcap","ceradrcg_hcap",
"animfl_hcap","bmim_hcap","bmdel_hcap","lmim_hcap","lmdel_hcap","lmrcg_hcap",
"conprxim_hcap","conprxdel_hcap","dsmt_hcap","numser_hcap","ravens_hcap",
"trailsa_hcap","trailsb_hcap","spatial_hcap")
### Recode MMSE
# summary(hcog[,mmse_nms])
hcog$imrec_mmse <- missCode(hcog,"imrec_mmse",0,3)
hcog$delrec_mmse <- missCode(hcog,"delrec_mmse",0,3)
# table(hcog$ornttime_mmse)
# table(hcog$orntplace_mmse)
# table(hcog$imrec_mmse)
# table(hcog$delrec_mmse)
# table(hcog$spellbck_mmse)
# table(hcog$naming_mmse)
# table(hcog$phrase_mmse)
# table(hcog$commfoll_mmse)
# table(hcog$writesent_mmse)
# table(hcog$draw_mmse)
# table(hcog$comm3step_mmse)
hcog <- hcog %>% mutate(
ornttime_mmser = ornttime_mmse,
orntplace_mmser = orntplace_mmse,
imrec_mmser = case_when(
imrec_mmse %in% c(0,1) ~ 0,
TRUE ~ imrec_mmse - 1
),
delrec_mmser = delrec_mmse,
spellbck_mmser = spellbck_mmse,
naming_mmser = naming_mmse,
phrase_mmser = phrase_mmse,
commfoll_mmser = commfoll_mmse,
writesent_mmser = writesent_mmse,
draw_mmser = draw_mmse,
comm3step_mmser = comm3step_mmse
)
# table(hcog$imrec_mmse,hcog$imrec_mmser)
### End recode MMSE
### Recode TICS
# summary(hcog[,tics_nms])
#
# table(hcog$wlimrc_tics)
# table(hcog$wldlrc_tics)
# table(hcog$cntbck_tics)
# table(hcog$count_tics)
# table(hcog$animfl_tics)
# table(hcog$numser_tics)
# table(hcog$serial7_tics)
# table(hcog$ornttime_tics)
# table(hcog$naming_tics)
hcog <- hcog %>% mutate(
wlimrc_ticsr = case_when(
wlimrc_tics %in% c(9,10) ~ 9,
TRUE ~ wlimrc_tics
),
wldlrc_ticsr = case_when(
wldlrc_tics %in% c(9,10) ~ 9,
TRUE ~ wldlrc_tics
),
cntbck_ticsr = case_when(
cntbck_tics %in% c(0,1) ~ 0,
TRUE ~ cntbck_tics - 1
),
# count_ticsr = case_when(
# count_tics %in% c(0,1,2) ~ 0,
# TRUE ~ count_tics - 2
# ),
ornttime_ticsr = case_when(
ornttime_tics %in% c(0,1) ~ 0,
TRUE ~ ornttime_tics - 1
),
naming_ticsr = case_when(
naming_tics %in% c(0,1,2) ~ 0,
TRUE ~ naming_tics - 2
),
numser_ticsr = numser_tics,
serial7_ticsr = serial7_tics
)
# table(hcog$wlimrc_tics,hcog$wlimrc_ticsr)
# table(hcog$wldlrc_tics,hcog$wldlrc_ticsr)
# table(hcog$count_tics,hcog$count_ticsr)
# table(hcog$ornttime_tics,hcog$ornttime_ticsr)
# table(hcog$naming_tics,hcog$naming_ticsr)
hcog <- recodeOrdinal(df = hcog, varlist_orig = "animfl_tics", varlist_tr = "animfl_ticsr")
# table(hcog$animfl_tics,hcog$animfl_ticsr)
### End recode TICS
### Recode HCAP
# summary(hcog[,hcap_nms])
hcog$trailsa_hcap <- missCode(hcog,"trailsa_hcap",0,300)
hcog$trailsb_hcap <- missCode(hcog,"trailsb_hcap",0,300)
hcog$trailsa_hcapi <- -hcog$trailsa_hcap
hcog$trailsb_hcapi <- -hcog$trailsb_hcap
# summary(hcog[,c("trailsa_hcapi","trailsb_hcapi")])
#
# table(hcog$naming_hcap)
# table(hcog$ceradtot_hcap)
# table(hcog$ceraddel_hcap)
# table(hcog$ceradrcg_hcap)
# table(hcog$animfl_hcap)
# table(hcog$bmim_hcap)
# table(hcog$bmdel_hcap)
# table(hcog$lmim_hcap)
# table(hcog$lmdel_hcap)
# table(hcog$lmrcg_hcap)
# table(hcog$conprxim_hcap)
# table(hcog$conprxdel_hcap)
# table(hcog$dsmt_hcap)
# table(hcog$numser_hcap)
# table(hcog$ravens_hcap)
# table(hcog$trailsa_hcap)
# table(hcog$trailsb_hcap)
# table(hcog$spatial_hcap)
hcog <- hcog %>% mutate(
naming_hcapr = case_when(
naming_hcap %in% c(0,1,2,3) ~ 1,
TRUE ~ naming_hcap - 2
),
spatial_hcapr = spatial_hcap + 1
)
varlist_orig <- c("ceradtot_hcap","ceraddel_hcap","ceradrcg_hcap",
"animfl_hcap","bmim_hcap","bmdel_hcap","lmim_hcap","lmdel_hcap","lmrcg_hcap",
"conprxim_hcap","conprxdel_hcap","dsmt_hcap","numser_hcap","ravens_hcap",
"trailsa_hcapi","trailsb_hcapi")
varlist_tr <- paste0(varlist_orig,"r")
hcog <- recodeOrdinal(hcog, varlist_orig = varlist_orig, varlist_tr = varlist_tr)
# table(hcog$naming_hcapr)
# table(hcog$ceradtot_hcapr)
# table(hcog$ceraddel_hcapr)
# table(hcog$ceradrcg_hcapr)
# table(hcog$animfl_hcapr)
# table(hcog$bmim_hcapr)
# table(hcog$bmdel_hcapr)
# table(hcog$lmim_hcapr)
# table(hcog$lmdel_hcapr)
# table(hcog$lmrcg_hcapr)
# table(hcog$conprxim_hcapr)
# table(hcog$conprxdel_hcapr)
# table(hcog$dsmt_hcapr)
# table(hcog$numser_hcapr)
# table(hcog$ravens_hcapr)
# table(hcog$trailsa_hcapir)
# table(hcog$trailsb_hcapir)
# table(hcog$spatial_hcapr)
### End recode HCAP
saveRDS(hcog, file = "~/Psychometrics Conference/2024/Simulation WG/Data/hrs_hcap_co_calibration_analytic.rds")
write.csv(hcog, file = "~/Psychometrics Conference/2024/Simulation WG/Data/hrs_hcap_co_calibration_analytic.csv", row.names=FALSE)
# --------------------------------- End Data --------------------------------
Multidimensional Item Response Theory Models
# ******************************* mirt models ********************************
hcog <- readRDS("~/Psychometrics Conference/2024/Simulation WG/Data/hrs_hcap_co_calibration_analytic.rds")
### Item names for analysis
mmse_nmsr <- c("ornttime_mmser","orntplace_mmser","imrec_mmser","delrec_mmser",
"spellbck_mmser","naming_mmser","phrase_mmser","commfoll_mmser","writesent_mmser",
"draw_mmser","comm3step_mmser")
tics_nmsr <- c("wlimrc_ticsr","wldlrc_ticsr","cntbck_ticsr","animfl_ticsr",
"numser_ticsr","serial7_ticsr","ornttime_ticsr","naming_ticsr")
hcap_nmsr <- c("naming_hcapr","ceradtot_hcapr","ceraddel_hcapr","ceradrcg_hcapr",
"animfl_hcapr","bmim_hcapr","bmdel_hcapr","lmim_hcapr","lmdel_hcapr","lmrcg_hcapr",
"conprxim_hcapr","conprxdel_hcapr","dsmt_hcapr","numser_hcapr","ravens_hcapr",
"trailsa_hcapir","trailsb_hcapir","spatial_hcapr")
hcog1 <- hcog[,c(mmse_nmsr,tics_nmsr,hcap_nmsr)]
names(hcog1)
### Unidimensional model
model_hrs_hcap_mmse <- "cog = 1-37"
res_1f <- mirt(hcog1, model_hrs_hcap_mmse,
method="MHRM",
technical=list(NCYCLES=1000))
summary(res_1f)
fit_1f <- M2(res_1f,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE)
### Multidimensional model - shared methods
model_hrs_hcap_mmse_m1 <- "cog = 1-37
recmmse = 3-4
wltics = 12-13
cerhcap = 21-23
bmhcap = 25-26
lmhcap = 27-29
cnpxhcap = 30-31
trhcap = 35-36"
res_md_1 <- mirt(hcog1, model_hrs_hcap_mmse_m1,
method="MHRM",
technical=list(NCYCLES=1000))
summary(res_md_1)
coef(res_md_1)
# par <- mod2values(res_md_1)
fit_md_1 <- M2(res_md_1,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE,
QMC = TRUE)
# constrained loadings of 2-indicator factors
res_md_1a <- mirt(hcog1, model_hrs_hcap_mmse_m1,
method="MHRM",
constrain = list(c(28,38),c(128,145),c(310,327),c(396,412),c(479,495)),
technical=list(NCYCLES=1000))
summary(res_md_1a)
coef(res_md_1a)
fit_md_1a <- M2(res_md_1a,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE,
QMC = TRUE)
### Multidimensional model - shared methods and content
model_hrs_hcap_mmse_m2 <- "cog = 1-37
recmmse = 3-4
wltics = 12-13
cerhcap = 21-23
bmhcap = 25-26
lmhcap = 27-29
cnpxhcap = 30-31
trhcap = 35-36
ornt = 1,18
nmng = 6,19,20"
res_md_2 <- mirt(hcog1, model_hrs_hcap_mmse_m2,
method="MHRM",
technical=list(NCYCLES=1000))
summary(res_md_2)
coef(res_md_2)
fit_md_2 <- M2(res_md_2,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE,
QMC = TRUE)
# constrained loadings of 2-indicator factors
res_md_2a <- mirt(hcog1, model_hrs_hcap_mmse_m2,
method="MHRM",
constrain = list(c(32,44),c(150,169),c(358,377),c(454,472),c(547,565),c(9,254)),
technical=list(NCYCLES=1000))
# par <- mod2values(res_md_2)
summary(res_md_2a)
coef(res_md_2a)
fit_md_2a <- M2(res_md_2a,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE,
QMC = TRUE)
save(res_1f,res_md_1,res_md_1a,res_md_2,res_md_2a,fit_1f,fit_md_1,fit_md_1a,fit_md_2,fit_md_2a,
file = "~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/co_calibration_results.Rdata")
# rm(list = ls()[grepl("^fit_",ls())])
# ------------------------------end mirt models ------------------------------
Model Fit Comparisons
load("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/co_calibration_results.Rdata")
fit_1f <- fit_1f %>% mutate(model = "unidimensional") %>%
relocate(model)
fit_md_1a <- fit_md_1a %>% mutate(model = "multidimensional 1") %>%
relocate(model)
fit_md_2a <- fit_md_2a %>% mutate(model = "multidimensional 2") %>%
relocate(model)
fit_summ <- bind_rows(fit_1f, fit_md_1a, fit_md_2a)
pandoc.table(fit_summ, row.names = FALSE, caption = "Model fit comparisons")
Model fit comparisons (continued below)
| unidimensional |
4100 |
629 |
0 |
0.06927 |
0.06723 |
0.07127 |
0.1571 |
| multidimensional 1 |
1993 |
618 |
0 |
0.04399 |
0.04182 |
0.04614 |
0.1493 |
| multidimensional 2 |
1991 |
614 |
0 |
0.04416 |
0.04199 |
0.04632 |
0.1488 |
| 0.8226 |
0.8325 |
| 0.9285 |
0.9336 |
| 0.9279 |
0.9335 |
Test Information Curves
extractMIRTParm <- function(model,max_thresh = 9) {
require(mirt)
parm <- mod2values(model)
parmw <- parm %>% dplyr::select(item, name, value) %>%
pivot_wider(id_cols = item, names_from = name, values_from = value) %>%
filter(!item == "GROUP") %>% mutate(
d1 = case_when(
is.na(d1) ~ d,
TRUE ~ d1
)
)
a <- as.matrix(parmw %>% dplyr::select(a1))
d <- as.matrix(parmw %>% dplyr::select(paste0("d",1:max_thresh)))
return(list(a,d))
}
hcog <- readRDS("~/Psychometrics Conference/2024/Simulation WG/Data/hrs_hcap_co_calibration_analytic.rds")
### Item names for analysis
mmse_nmsr <- c("ornttime_mmser","orntplace_mmser","imrec_mmser","delrec_mmser",
"spellbck_mmser","naming_mmser","phrase_mmser","commfoll_mmser","writesent_mmser",
"draw_mmser","comm3step_mmser")
tics_nmsr <- c("wlimrc_ticsr","wldlrc_ticsr","cntbck_ticsr","animfl_ticsr",
"numser_ticsr","serial7_ticsr","ornttime_ticsr","naming_ticsr")
hcap_nmsr <- c("naming_hcapr","ceradtot_hcapr","ceraddel_hcapr","ceradrcg_hcapr",
"animfl_hcapr","bmim_hcapr","bmdel_hcapr","lmim_hcapr","lmdel_hcapr","lmrcg_hcapr",
"conprxim_hcapr","conprxdel_hcapr","dsmt_hcapr","numser_hcapr","ravens_hcapr",
"trailsa_hcapir","trailsb_hcapir","spatial_hcapr")
hcog1 <- hcog[,c(mmse_nmsr,tics_nmsr,hcap_nmsr)]
mirt_mmse <- "mmse = 1-11"
mirt_tics3 <- "tics = 1-6"
mirt_tics10 <- "tics = 1-7"
mirt_tics13 <- "tics = 1-8"
mirt_tics16 <- "tics = 1-10"
mirt_hcap <- "hcap = 1-18"
mirt_all <- "cog = 1-37"
true_sim <- readRDS("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Data/simulated_longitudinal_true_cognition.rds")
item_labels <- names(hcog1)
a <- extractMIRTParm((res_md_1a))[[1]]
d <- extractMIRTParm((res_md_1a))[[2]]
df <- data.frame(simdata(a = a, d = d, N = 1000, itemtype = 'graded'))
# cat("Test Information Curve - MMSE")
mod <- mirt(data = df[,1:11], model = mirt_mmse)
info_mmse <- plot(mod, type = 'info', ylim = c(0,40),
main = "Test Information Curve - MMSE")
# cat("Test Information Curve - TICS3")
mod <- mirt(data = df[,c(12:14,17:19)], model = mirt_tics3)
info_tics3 <- plot(mod, type = 'info', ylim = c(0,40),
main = "Test Information Curve - TICS3")
# cat("Test Information Curve - TICS10")
mod <- mirt(data = df[,c(12:14,16:19)], model = mirt_tics10)
info_tics10 <- plot(mod, type = 'info', ylim = c(0,40),
main = "Test Information Curve - TICS10")
# cat("Test Information Curve - TICS13")
mod <- mirt(data = df[,12:19], model = mirt_tics13)
info_tics13 <- plot(mod, type = 'info', ylim = c(0,40),
main = "Test Information Curve - TICS13")
# cat("Test Information Curve - TICS16")
mod <- mirt(data = df[,c(12:19,35:36)], model = mirt_tics16)
info_tics16 <- plot(mod, type = 'info', ylim = c(0,40),
main = "Test Information Curve - TICS16")
# cat("Test Information Curve - HCAP")
mod <- mirt(data = df[,20:37], model = mirt_hcap)
info_hcap <- plot(mod, type = 'info', ylim = c(0,40),
main = "Test Information Curve - HCAP")
# cat("Test Information Curve - MMSE + TICS + HCAP")
mod <- mirt(data = df, model = mirt_all)
info_all <- plot(mod, type = 'info', ylim = c(0,40),
main = "Test Information Curve - MMSE + TICS + HCAP")
Simulation of Longitudinal Cognitive Change
Simulation Functions
require(lme4)
require(mirt)
require(tidyverse)
require(ggplot2)
# show size of environment objects
# sort( sapply(ls(),function(x){object.size(get(x))}))
# rm(list = ls()[grepl("^re",ls())])
# extractMIRTParm is a function that extracts item parameters from a mirt
# estimate multidimensional (uncorrelated dimensions) model and converts
# discrimination (a) parameters to a vector of item parameters and
# threshhold (d) parameters to a matrix of d values.
# Parameters
# model - estimated mirt model object
# max_thresh - maximum number of threshold paramters across items
# Value list with 2 elements
# a - vector of a parameters, 1 for each item
# d - matrix of d parameters, rows correspons to items, columns to threshholds
extractMIRTParm <- function(model,max_thresh = 9) {
require(mirt)
parm <- mod2values(model)
parmw <- parm %>% dplyr::select(item, name, value) %>%
pivot_wider(id_cols = item, names_from = name, values_from = value) %>%
filter(!item == "GROUP") %>% mutate(
d1 = case_when(
is.na(d1) ~ d,
TRUE ~ d1
)
)
a <- as.matrix(parmw %>% dplyr::select(a1))
d <- as.matrix(parmw %>% dplyr::select(paste0("d",1:max_thresh)))
return(list(a,d))
}
# simTraj is a function that 1) estimates a mirt model using a simulated
# item response dataset, 2) calculates factor scores based on the estimated
# model, 3) estimates linear mixed effects models with true cognition and
# simulated cognition (factor scores from simulated item responses) as dependent
# variables, time in study as a fixed effect variable, and person and
# person-by-time random effects, and 4) returns person specific estimates
# (random effects) of cognitive components slopes and intercepts.
# Parameters
# data - data frame that contains true cognition value, simulated item responses,
# id variable, and time in study variable
# sim_cog_var - label for variable within dataframe (data) for which
# trajectories are being simulated
# frml - formula specification for (lmer) longitudinal mixed effects model
# iteration = iteration number passed from simulateTrajectories
# Value - list with 3 elements:
# data - analytic dataframe (data) with longitudinal model predicted cognition
# values for each assessment
# re - data frame with estimated random effects from longitudinal model.
# Includes intercept and slope random effects for true cognition (_true)
# and simulated cognition factor scores (_sim)
# fe - data frame with estimated fixed effects from longitudinal model.
# simTraj <- function(data = df, item_labels = item_labels, mirt_mod = mirt_all,
# iteration = iteration) {
simTraj <- function(data = df, sim_cog_var = "sim_cog_all", frml = frml,
iteration = iteration) {
require(tidyverse)
require(lme4)
# # estimate mirt model using simulated responses and generate factor scores
re_empty <- data.frame(id = NA, int_true = NA, slope_true = NA,
int_true_se = NA,slope_true_se = NA, int_sim = NA, slope_sim = NA,
int_sim_se = NA, slope_sim_se = NA)
#
# tryCatch({
# mod <- mirt(data[,item_labels], mirt_mod)
# data$sim_cog <- fscores(mod)
# data <- data %>% relocate(sim_cog, .after = true_cog)
# }, error=function(e){return(re_empty)})
# data$itertion <- iteration
# frml <- "true_cog ~ time + agebl_75 + slope_true_qrtl +
# time:agebl_75 + time:slope_true_qrtl + (1 + time | id)"
# frml <- "true_cog ~ time + agebl_75 +
# time:agebl_75 + (1 | id)"
#
tryCatch({
# mixed effects longitudinal model - true cognition
# res_long_true <- lmer(true_cog ~ time + (1 + time | id), data = data)
res_long_true <- lmer(formula = frml, data = data)
# summary(res_long_true)
# str(coef(res_long_true))
data$cog_true_pred <- predict(res_long_true, newdata = data, type = "response")
re_true <- data.frame(ranef(res_long_true))
re_true <- re_true %>%
pivot_wider(id_cols = grp, names_from = term,values_from = c(condval,condsd))
names(re_true) <- c("id","int_true","slope_true","int_true_se","slope_true_se")
re_true$id <- as.integer(levels(re_true$id))[re_true$id]
fe_true <- data.frame(t(lme4::fixef(res_long_true)))
names(fe_true) <- paste0(names(fe_true),"_true")
data$sim_cog <- data[,sim_cog_var]
#res_long_sim <- lmer(sim_cog ~ time + (1 + time | id), data = data)
frml <- sub("true_cog","sim_cog",frml)
res_long_sim <- lmer(formula = frml, data = data)
# summary(res_long_sim)
data$cog_sim_pred <- predict(res_long_sim, newdata = data, type = "response")
data <- data %>% dplyr::select(-sim_cog)
re_sim <- data.frame(ranef(res_long_sim))
re_sim <- re_sim %>%
pivot_wider(id_cols = grp, names_from = term,values_from = c(condval,condsd))
names(re_sim) <- c("id","int_sim","slope_sim","int_sim_se","slope_sim_se")
re_sim$id <- as.integer(levels(re_sim$id))[re_sim$id]
fe_sim <- data.frame(t(lme4::fixef(res_long_sim)))
names(fe_sim) <- paste0(names(fe_sim),"_sim")
re <- re_true %>% left_join(re_sim, by="id")
fe <- bind_cols(fe_true,fe_sim)
fe$iteration <- iteration
}, error=function(e){re <- re_empty})
re$iteration <- iteration
return(list(data,re,fe))
}
# simulateTrajectories is a function that 1) simulates item responses for
# true cognition vectors for different components of the HRS-HCAP cognitive test
# battery (MMSE, HRS TICS, HCAP, MMSE+TICS+HCAP), and 2) generates simulated observed (factor)
# scores for each of the cognitive components (simTraj()), 3) estimates linear mixed effects models
# with true cognition and simulated cognition (factor scores from simulated item responses)
# as dependent variables, time in study as a fixed effect variable, and person and
# person-by-time random effects (simTraj()), 4) collates person specific estimates (random effects)
# of cognitive components slopes and intercepts.
#
# Parameters
# true_sim - data frame with true cognition values, rows correspond to assessments,
# columns correspond to simulation iterations.
# a_par - vector/data frame containing a (discrimination) item parameters
# from mirt co-calibration model
# d_par - vector/data frame containing d (threshhold) item parameters
# from mirt co-calibration model
# item_labels - item labels
# iter_group - number of groups of (niter) iterations
# niter - number of simulation iterations within iteration groups
# out_dir - path to directory for output of simulation results
# frml - formula specification for (lmer) longitudinal mixed effects model
#
# Value - Saves lists of results to out_dir - there is a file for each iter_group
# that contains niter elements:
# ds - list of datasets for each niter simulations. Each dataset contains
# the true cognition value (true_cog), the simulated item responses conditional
# on true_cog, simulated item response data and cognition factor scores,
# and longitudinal mixed effect model predicted cognition values for each simulated
# assessment. Measures include MMSE, TICS, HCAP, and MMSE+TICS+HCAP.
# re_all - estimated random effects for the full HRS-HCAP cognitive test battery
# (MMSE+TICS+HCAP). Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim)
# re_mmse - estimated random effects for the MMSE component of the HRS-HCAP
# cognitive test battery. Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim)
# re_tics - estimated random effects for the TICS component of the HRS-HCAP
# cognitive test battery. Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim)
# re_hcap - estimated random effects for the HCAP component of the HRS-HCAP
# cognitive test battery. Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim)
# fe_all - estimated fixed effects for mixed effect model of cognition measured
# by full HRS-HCAP cognitive test battery (MMSE+TICS+HCAP).
# Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim).
# fe_mmse - estimated fixed effects for mixed effect model of cognition measured
# by MMSE.
# Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim).
# fe_tics - estimated fixed effects for mixed effect model of cognition measured
# by TICS.
# Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim).
# fe_hcap - estimated fixed effects for mixed effect model of cognition measured
# by HCAP.
# Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim).
simulateTrajectories <- function(true_sim, a_par, d_par, item_labels,
niter = 20, iter_group = 5, out_dir = "Analysis/Simulation Results/",
frml = "true_cog ~ time + (1 | id)") {
require(mirt)
require(tidyverse)
for (itgrp in 1:iter_group) {
ds <- list()
re_all <- list()
re_mmse <- list()
re_tics <- list()
re_hcap <- list()
fe_all <- list()
fe_mmse <- list()
fe_tics <- list()
fe_hcap <- list()
set.seed(092724)
for (iter in 1:niter) {
cat(paste0("Group - ",itgrp, ", Iteration - ",iter,"\n"))
# simulate item responses
iteration <- ((itgrp - 1) * niter) + iter
df <- data.frame(simdata(a = a_par, d = d_par, Theta = true_sim[,paste0("vm",iteration)],
itemtype = 'graded'))
names(df) <- item_labels
df$id <- true_sim$id
df$time <- true_sim$time
df$agebl_75 <- true_sim$agebl_75
df$true_cog <- true_sim[,paste0("vm",iteration)]
df$iteration <- iteration
df <- df %>% relocate(c(id,iteration,time,true_cog))
# true random effect intercept and slope - from calculate_re_true.R
true_re <- readRDS("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Data/true_random_effects.rds")
df <- df %>% left_join((true_re %>% dplyr::select(id,int_true_qrtl,
slope_true_qrtl)),by="id")
flag_low_response <- FALSE
for (itnm in item_labels) {
if (length(table(df[,itnm])) < 2) {
flag_low_response <- TRUE
}
}
if (flag_low_response == TRUE) {
iter = iter - 1
next
}
# mirt_models for cognition measures
mirt_mmse <- "mmse = 1-11"
mirt_tics <- "tics = 1-8"
mirt_hcap <- "hcap = 1-18"
mirt_all <- "cog = 1-37"
### estimate mirt model using simulated responses and generate factor scores
# create empty re dataframe to return if error in generation
re_empty <- data.frame(id = NA, int_true = NA, slope_true = NA,
int_true_se = NA,slope_true_se = NA, int_sim = NA, slope_sim = NA,
int_sim_se = NA, slope_sim_se = NA)
tryCatch({
mod <- mirt(df[,item_labels], mirt_all)
df$sim_cog_all <- fscores(mod)
mod <- mirt(df[,item_labels[1:11]], mirt_mmse)
df$sim_cog_mmse <- fscores(mod)
mod <- mirt(df[,item_labels[12:19]], mirt_tics)
df$sim_cog_tics <- fscores(mod)
mod <- mirt(df[,item_labels[20:37]], mirt_hcap)
df$sim_cog_hcap <- fscores(mod)
df <- df %>%
relocate(c(sim_cog_all,sim_cog_mmse,sim_cog_tics,sim_cog_hcap),
.after = true_cog)
}, error=function(e){return(re_empty)})
# rescale factor scores to equate metric to true cognition
df <- df %>% mutate(
sim_cog_mmse_rs = (((sim_cog_mmse - mean(sim_cog_mmse)) /
sd(sim_cog_mmse)) * sd(true_cog)) + mean(true_cog),
sim_cog_tics_rs = (((sim_cog_tics - mean(sim_cog_tics)) /
sd(sim_cog_tics)) * sd(true_cog)) + mean(true_cog),
sim_cog_hcap_rs = (((sim_cog_hcap - mean(sim_cog_hcap)) /
sd(sim_cog_hcap)) * sd(true_cog)) + mean(true_cog),
sim_cog_all_rs = (((sim_cog_all - mean(sim_cog_all)) /
sd(sim_cog_all)) * sd(true_cog)) + mean(true_cog),
)
### calculate simulated random effects
# frml <- as.formula(frml)
# frml <- as.formula("true_cog ~ time + (1 | id)")
# frml <- "true_cog ~ time + agebl_75 + slope_true_qrtl +
# time:agebl_75 + time:slope_true_qrtl + (1 + time | id)"
# MMSE+TICS+HCAP
res <- simTraj(data = df, sim_cog_var = "sim_cog_all_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_all = cog_true_pred,
cog_sim_pred_all = cog_sim_pred)
re <- res[[2]]
re$label <- "MMSE+TICS+HCAP"
re_all[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "MMSE+TICS+HCAP"
fe_all[[paste0("iteration-",iteration)]] <- fe
# MMSE
res <- simTraj(data = df, sim_cog_var = "sim_cog_mmse_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_mmse = cog_true_pred,
cog_sim_pred_mmse = cog_sim_pred)
re <- res[[2]]
re$label <- "MMSE"
re_mmse[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "MMSE"
fe_mmse[[paste0("iteration-",iteration)]] <- fe
# TICS
res <- simTraj(data = df, sim_cog_var = "sim_cog_tics_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_tics = cog_true_pred,
cog_sim_pred_tics = cog_sim_pred)
re <- res[[2]]
re$label <- "TICS"
re_tics[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "TICS"
fe_tics[[paste0("iteration-",iteration)]] <- fe
# HCAP
res <- simTraj(data = df, sim_cog_var = "sim_cog_hcap_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_hcap = cog_true_pred,
cog_sim_pred_hcap = cog_sim_pred)
re <- res[[2]]
re$label <- "HCAP"
re_hcap[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "HCAP"
fe_hcap[[paste0("iteration-",iteration)]] <- fe
ds[[paste0("iteration-",iteration)]] <- df
} # end for iter
saveRDS(ds, file = paste0(out_dir,"ds_iteration_group_", itgrp, ".rds"))
saveRDS(re_mmse, file = paste0(out_dir,"re_mmse_iteration_group_", itgrp, ".rds"))
saveRDS(re_tics, file = paste0(out_dir,"re_tics_iteration_group_", itgrp, ".rds"))
saveRDS(re_hcap, file = paste0(out_dir,"re_hcap_iteration_group_", itgrp, ".rds"))
saveRDS(re_all, file = paste0(out_dir,"re_all_iteration_group_", itgrp, ".rds"))
saveRDS(fe_mmse, file = paste0(out_dir,"fe_mmse_iteration_group_", itgrp, ".rds"))
saveRDS(fe_tics, file = paste0(out_dir,"fe_tics_iteration_group_", itgrp, ".rds"))
saveRDS(fe_hcap, file = paste0(out_dir,"fe_hcap_iteration_group_", itgrp, ".rds"))
saveRDS(fe_all, file = paste0(out_dir,"fe_all_iteration_group_", itgrp, ".rds"))
} # end for itgrp
for (itgrp in 1:iter_group) {
if (itgrp == 1) {
re_mmse <- readRDS(paste0(out_dir,"re_mmse_iteration_group_", itgrp, ".rds"))
re_mmse <- re_mmse %>% bind_rows()
re_tics <- readRDS(paste0(out_dir,"re_tics_iteration_group_", itgrp, ".rds"))
re_tics <- re_tics %>% bind_rows()
re_hcap <- readRDS(paste0(out_dir,"re_hcap_iteration_group_", itgrp, ".rds"))
re_hcap <- re_hcap %>% bind_rows()
re_all <- readRDS(paste0(out_dir,"re_all_iteration_group_", itgrp, ".rds"))
re_all <- re_all %>% bind_rows()
fe_mmse <- readRDS(paste0(out_dir,"fe_mmse_iteration_group_", itgrp, ".rds"))
fe_mmse <- fe_mmse %>% bind_rows()
fe_tics <- readRDS(paste0(out_dir,"fe_tics_iteration_group_", itgrp, ".rds"))
fe_tics <- fe_tics %>% bind_rows()
fe_hcap <- readRDS(paste0(out_dir,"fe_hcap_iteration_group_", itgrp, ".rds"))
fe_hcap <- fe_hcap %>% bind_rows()
fe_all <- readRDS(paste0(out_dir,"fe_all_iteration_group_", itgrp, ".rds"))
fe_all <- fe_all %>% bind_rows()
} else {
re_mmse_temp <- readRDS(paste0(out_dir,"re_mmse_iteration_group_", itgrp, ".rds"))
re_mmse_temp <- re_mmse_temp %>% bind_rows()
re_mmse <- bind_rows(re_mmse, re_mmse_temp)
re_tics_temp <- readRDS(paste0(out_dir,"re_tics_iteration_group_", itgrp, ".rds"))
re_tics_temp <- re_tics_temp %>% bind_rows()
re_tics <- bind_rows(re_tics, re_tics_temp)
re_hcap_temp <- readRDS(paste0(out_dir,"re_hcap_iteration_group_", itgrp, ".rds"))
re_hcap_temp <- re_hcap_temp %>% bind_rows()
re_hcap <- bind_rows(re_hcap, re_hcap_temp)
re_all_temp <- readRDS(paste0(out_dir,"re_all_iteration_group_", itgrp, ".rds"))
re_all_temp <- re_all_temp %>% bind_rows()
re_all <- bind_rows(re_all, re_all_temp)
fe_mmse_temp <- readRDS(paste0(out_dir,"fe_mmse_iteration_group_", itgrp, ".rds"))
fe_mmse_temp <- fe_mmse_temp %>% bind_rows()
fe_mmse <- bind_rows(fe_mmse, fe_mmse_temp)
fe_tics_temp <- readRDS(paste0(out_dir,"fe_tics_iteration_group_", itgrp, ".rds"))
fe_tics_temp <- fe_tics_temp %>% bind_rows()
fe_tics <- bind_rows(fe_tics, fe_tics_temp)
fe_hcap_temp <- readRDS(paste0(out_dir,"fe_hcap_iteration_group_", itgrp, ".rds"))
fe_hcap_temp <- fe_hcap_temp %>% bind_rows()
fe_hcap <- bind_rows(fe_hcap, fe_hcap_temp)
fe_all_temp <- readRDS(paste0(out_dir,"fe_all_iteration_group_", itgrp, ".rds"))
fe_all_temp <- fe_all_temp %>% bind_rows()
fe_all <- bind_rows(fe_all, fe_all_temp)
}
}
saveRDS(re_mmse, file = paste0(out_dir,"re_mmse", ".rds"))
saveRDS(re_tics, file = paste0(out_dir,"re_tics", ".rds"))
saveRDS(re_hcap, file = paste0(out_dir,"re_hcap", ".rds"))
saveRDS(re_all, file = paste0(out_dir,"re_all", ".rds"))
saveRDS(fe_mmse, file = paste0(out_dir,"fe_mmse", ".rds"))
saveRDS(fe_tics, file = paste0(out_dir,"fe_tics", ".rds"))
saveRDS(fe_hcap, file = paste0(out_dir,"fe_hcap", ".rds"))
saveRDS(fe_all, file = paste0(out_dir,"fe_all", ".rds"))
# return(list("ds" = ds, "re_all" = re_all, "re_mmse" = re_mmse,
# "re_tics" = re_tics, "re_hcap" = re_hcap))
}
simulateTrajectories2 <- function(item_labels, niter = 20, iter_group = 5,
in_dir = "Analysis/Simulation Results/",
out_dir = "Analysis/Simulation Results/",
frml = "true_cog ~ time + (1 | id)") {
require(mirt)
require(tidyverse)
for (itgrp in 1:iter_group) {
ds <- list()
re_tics3 <- list()
re_tics10 <- list()
re_tics13 <- list()
re_tics16 <- list()
fe_tics3 <- list()
fe_tics10 <- list()
fe_tics13 <- list()
fe_tics16 <- list()
set.seed(092724)
for (iter in 1:niter) {
cat(paste0("Group - ",itgrp, ", Iteration - ",iter,"\n"))
# load simulated item responses
iteration <- ((itgrp - 1) * niter) + iter
df <- readRDS(paste0(in_dir,"ds_iteration_group_",itgrp,".rds"))[[iter]]
df <- df %>% dplyr::select(c(id:true_cog,agebl_75:slope_true_qrtl,any_of(item_labels)))
# ds_iteration_group_36.rds
#
# df <- data.frame(simdata(a = a_par, d = d_par, Theta = true_sim[,paste0("vm",iteration)],
# itemtype = 'graded'))
# names(df) <- item_labels
# df$id <- true_sim$id
# df$time <- true_sim$time
# df$agebl_75 <- true_sim$agebl_75
# df$true_cog <- true_sim[,paste0("vm",iteration)]
# df$iteration <- iteration
# df <- df %>% relocate(c(id,iteration,time,true_cog))
# # true random effect intercept and slope - from calculate_re_true.R
# true_re <- readRDS("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Data/true_random_effects.rds")
# df <- df %>% left_join((true_re %>% dplyr::select(id,int_true_qrtl,
# slope_true_qrtl)),by="id")
#
# flag_low_response <- FALSE
# for (itnm in item_labels) {
# if (length(table(df[,itnm])) < 2) {
# flag_low_response <- TRUE
# }
# }
# if (flag_low_response == TRUE) {
# iter = iter - 1
# next
# }
# mirt_models for cognition measures
# mirt_mmse <- "mmse = 1-11"
# mirt_tics <- "tics = 1-8"
# mirt_hcap <- "hcap = 1-18"
# mirt_all <- "cog = 1-37"
mirt_tics3 <- "tics3 = 1-6"
mirt_tics10 <- "tics10 = 1-7"
mirt_tics13 <- "tics13 = 1-8"
mirt_tics16 <- "tics16 = 1-10"
### estimate mirt model using simulated responses and generate factor scores
# create empty re dataframe to return if error in generation
re_empty <- data.frame(id = NA, int_true = NA, slope_true = NA,
int_true_se = NA,slope_true_se = NA, int_sim = NA, slope_sim = NA,
int_sim_se = NA, slope_sim_se = NA)
tryCatch({
mod <- mirt(df[,item_labels[c(1:3,6:8)]], mirt_tics3)
df$sim_cog_tics3 <- fscores(mod)
mod <- mirt(df[,item_labels[c(1:3,5:8)]], mirt_tics10)
df$sim_cog_tics10 <- fscores(mod)
mod <- mirt(df[,item_labels[1:8]], mirt_tics13)
df$sim_cog_tics13 <- fscores(mod)
mod <- mirt(df[,item_labels[1:10]], mirt_tics16)
df$sim_cog_tics16 <- fscores(mod)
df <- df %>%
relocate(c(sim_cog_tics3,sim_cog_tics10,sim_cog_tics13,sim_cog_tics16),
.after = true_cog)
}, error=function(e){return(re_empty)})
# rescale factor scores to equate metric to true cognition
df <- df %>% mutate(
sim_cog_tics3_rs = (((sim_cog_tics3 - mean(sim_cog_tics3)) /
sd(sim_cog_tics3)) * sd(true_cog)) + mean(true_cog),
sim_cog_tics10_rs = (((sim_cog_tics10 - mean(sim_cog_tics10)) /
sd(sim_cog_tics10)) * sd(true_cog)) + mean(true_cog),
sim_cog_tics13_rs = (((sim_cog_tics13 - mean(sim_cog_tics13)) /
sd(sim_cog_tics13)) * sd(true_cog)) + mean(true_cog),
sim_cog_tics16_rs = (((sim_cog_tics16 - mean(sim_cog_tics16)) /
sd(sim_cog_tics16)) * sd(true_cog)) + mean(true_cog),
)
### calculate simulated random effects
# frml <- as.formula(frml)
# frml <- as.formula("true_cog ~ time + (1 | id)")
# frml <- "true_cog ~ time + agebl_75 + slope_true_qrtl +
# time:agebl_75 + time:slope_true_qrtl + (1 + time | id)"
# TICS3
res <- simTraj(data = df, sim_cog_var = "sim_cog_tics3_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_tics3 = cog_true_pred,
cog_sim_pred_tics3 = cog_sim_pred)
re <- res[[2]]
re$label <- "TICS3"
re_tics3[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "TICS3"
fe_tics3[[paste0("iteration-",iteration)]] <- fe
# TICS10
res <- simTraj(data = df, sim_cog_var = "sim_cog_tics10_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_tics10 = cog_true_pred,
cog_sim_pred_tics10 = cog_sim_pred)
re <- res[[2]]
re$label <- "TICS10"
re_tics10[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "TICS10"
fe_tics10[[paste0("iteration-",iteration)]] <- fe
# TICS13
res <- simTraj(data = df, sim_cog_var = "sim_cog_tics13_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_tics13 = cog_true_pred,
cog_sim_pred_tics13 = cog_sim_pred)
re <- res[[2]]
re$label <- "TICS13"
re_tics13[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "TICS13"
fe_tics13[[paste0("iteration-",iteration)]] <- fe
# TICS16
res <- simTraj(data = df, sim_cog_var = "sim_cog_tics16_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_tics16 = cog_true_pred,
cog_sim_pred_tics16 = cog_sim_pred)
re <- res[[2]]
re$label <- "TICS16"
re_tics16[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "TICS16"
fe_tics16[[paste0("iteration-",iteration)]] <- fe
ds[[paste0("iteration-",iteration)]] <- df
} # end for iter
saveRDS(ds, file = paste0(out_dir,"ds_iteration_group_", itgrp, ".rds"))
saveRDS(re_tics3, file = paste0(out_dir,"re_tics3_iteration_group_", itgrp, ".rds"))
saveRDS(re_tics10, file = paste0(out_dir,"re_tics10_iteration_group_", itgrp, ".rds"))
saveRDS(re_tics13, file = paste0(out_dir,"re_tics13_iteration_group_", itgrp, ".rds"))
saveRDS(re_tics16, file = paste0(out_dir,"re_tics16_iteration_group_", itgrp, ".rds"))
saveRDS(fe_tics3, file = paste0(out_dir,"fe_tics3_iteration_group_", itgrp, ".rds"))
saveRDS(fe_tics10, file = paste0(out_dir,"fe_tics10_iteration_group_", itgrp, ".rds"))
saveRDS(fe_tics13, file = paste0(out_dir,"fe_tics13_iteration_group_", itgrp, ".rds"))
saveRDS(fe_tics16, file = paste0(out_dir,"fe_tics16_iteration_group_", itgrp, ".rds"))
} # end for itgrp
for (itgrp in 1:iter_group) {
if (itgrp == 1) {
re_tics3 <- readRDS(paste0(out_dir,"re_tics3_iteration_group_", itgrp, ".rds"))
re_tics3 <- re_tics3 %>% bind_rows()
re_tics10 <- readRDS(paste0(out_dir,"re_tics10_iteration_group_", itgrp, ".rds"))
re_tics10 <- re_tics10 %>% bind_rows()
re_tics13 <- readRDS(paste0(out_dir,"re_tics13_iteration_group_", itgrp, ".rds"))
re_tics13 <- re_tics13 %>% bind_rows()
re_tics16 <- readRDS(paste0(out_dir,"re_tics16_iteration_group_", itgrp, ".rds"))
re_tics16 <- re_tics16 %>% bind_rows()
fe_tics3 <- readRDS(paste0(out_dir,"fe_tics3_iteration_group_", itgrp, ".rds"))
fe_tics3 <- fe_tics3 %>% bind_rows()
fe_tics10 <- readRDS(paste0(out_dir,"fe_tics10_iteration_group_", itgrp, ".rds"))
fe_tics10 <- fe_tics10 %>% bind_rows()
fe_tics13 <- readRDS(paste0(out_dir,"fe_tics13_iteration_group_", itgrp, ".rds"))
fe_tics13 <- fe_tics13 %>% bind_rows()
fe_tics16 <- readRDS(paste0(out_dir,"fe_tics16_iteration_group_", itgrp, ".rds"))
fe_tics16 <- fe_tics16 %>% bind_rows()
} else {
re_tics3_temp <- readRDS(paste0(out_dir,"re_tics3_iteration_group_", itgrp, ".rds"))
re_tics3_temp <- re_tics3_temp %>% bind_rows()
re_tics3 <- bind_rows(re_tics3, re_tics3_temp)
re_tics10_temp <- readRDS(paste0(out_dir,"re_tics10_iteration_group_", itgrp, ".rds"))
re_tics10_temp <- re_tics10_temp %>% bind_rows()
re_tics10 <- bind_rows(re_tics10, re_tics10_temp)
re_tics13_temp <- readRDS(paste0(out_dir,"re_tics13_iteration_group_", itgrp, ".rds"))
re_tics13_temp <- re_tics13_temp %>% bind_rows()
re_tics13 <- bind_rows(re_tics13, re_tics13_temp)
re_tics16_temp <- readRDS(paste0(out_dir,"re_tics16_iteration_group_", itgrp, ".rds"))
re_tics16_temp <- re_tics16_temp %>% bind_rows()
re_tics16 <- bind_rows(re_tics16, re_tics16_temp)
fe_tics3_temp <- readRDS(paste0(out_dir,"fe_tics3_iteration_group_", itgrp, ".rds"))
fe_tics3_temp <- fe_tics3_temp %>% bind_rows()
fe_tics3 <- bind_rows(fe_tics3, fe_tics3_temp)
fe_tics10_temp <- readRDS(paste0(out_dir,"fe_tics10_iteration_group_", itgrp, ".rds"))
fe_tics10_temp <- fe_tics10_temp %>% bind_rows()
fe_tics10 <- bind_rows(fe_tics10, fe_tics10_temp)
fe_tics13_temp <- readRDS(paste0(out_dir,"fe_tics13_iteration_group_", itgrp, ".rds"))
fe_tics13_temp <- fe_tics13_temp %>% bind_rows()
fe_tics13 <- bind_rows(fe_tics13, fe_tics13_temp)
fe_tics16_temp <- readRDS(paste0(out_dir,"fe_tics16_iteration_group_", itgrp, ".rds"))
fe_tics16_temp <- fe_tics16_temp %>% bind_rows()
fe_tics16 <- bind_rows(fe_tics16, fe_tics16_temp)
}
}
saveRDS(re_tics3, file = paste0(out_dir,"re_tics3", ".rds"))
saveRDS(re_tics10, file = paste0(out_dir,"re_tics10", ".rds"))
saveRDS(re_tics13, file = paste0(out_dir,"re_tics13", ".rds"))
saveRDS(re_tics16, file = paste0(out_dir,"re_tics16", ".rds"))
saveRDS(fe_tics3, file = paste0(out_dir,"fe_tics3", ".rds"))
saveRDS(fe_tics10, file = paste0(out_dir,"fe_tics10", ".rds"))
saveRDS(fe_tics13, file = paste0(out_dir,"fe_tics13", ".rds"))
saveRDS(fe_tics16, file = paste0(out_dir,"fe_tics16", ".rds"))
# return(list("ds" = ds, "re_all" = re_all, "re_mmse" = re_mmse,
# "re_tics" = re_tics, "re_hcap" = re_hcap))
}
# dat_cols <- c("id","time","age","agebl_75","slope_true_qrtl","int_true_qrtl")
# res_cols <- c("cog_true_pred_all","cog_sim_pred_all","cog_sim_pred_mmse",
# "cog_sim_pred_tics","cog_sim_pred_hcap")
mergeSimResults <-function(file_name = "ds_iteration_group_1.rds",
out_dir = out_dir,
dat_cols = dat_cols,
res_cols = res_cols,
iter_group = 50,
niter = 20) {
res <- (readRDS(paste0(out_dir,file_name))[[1]] %>%
dplyr::select(any_of(dat_cols)))
for (i in 1:iter_group) {
for (n in 1:niter) {
iteration <- ((i-1) * niter) + n
fl_nm <- sub("_1",paste0("_",i),file_name)
ds <- readRDS(paste0(out_dir,fl_nm))[[n]] %>%
dplyr::select(all_of(res_cols))
names(ds) <- paste0(res_cols,".",iteration)
res <- res %>% bind_cols(ds)
}
}
return(res)
}
Simulated Cognition Change versus True Cognition Change Plots
# sim20 <- readRDS("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/sim20.rds")
#
# re_all <- sim20[["re_all"]]
# re_mmse <- sim20[["re_mmse"]]
# re_tics <- sim20[["re_tics"]]
# re_hcap <- sim20[["re_hcap"]]
out_dir <- "~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/Simulation Results/10_04_24_unconditional/"
re_all <- readRDS(paste0(out_dir,"re_all.rds"))
re_mmse <- readRDS(paste0(out_dir,"re_mmse.rds"))
re_tics <- readRDS(paste0(out_dir,"re_tics.rds"))
re_hcap <- readRDS(paste0(out_dir,"re_hcap.rds"))
re_all_summ <- bind_rows(re_all)
re_mmse_summ <- bind_rows(re_mmse)
re_tics_summ <- bind_rows(re_tics)
re_hcap_summ <- bind_rows(re_hcap)
ggplot(data=re_mmse,aes(slope_true,slope_sim)) +
geom_point() +
geom_smooth() +
labs(caption = "True versus simulated cognitive change - MMSE")
ggplot(data=re_tics,aes(slope_true,slope_sim)) +
geom_point() +
geom_smooth() +
labs(caption = "True versus simulated cognitive change - TICS")
ggplot(data=re_hcap,aes(slope_true,slope_sim)) +
geom_point() +
geom_smooth() +
labs(caption = "True versus simulated cognitive change - HCAP")
ggplot(data=re_all,aes(slope_true,slope_sim)) +
geom_point() +
geom_smooth() +
labs(caption = "True versus simulated cognitive change - MMSE + TICS + HCAP")
Simulated Cognitive Trajectories versus True Cognition Trajectories
# Create model predicted trajectories
# # code to create merged r1 dataset
# dat_cols <- c("id","time","agebl_75","slope_true_qrtl","int_true_qrtl")
# res_cols <- c("cog_true_pred_all","cog_sim_pred_all","cog_sim_pred_mmse",
# "cog_sim_pred_tics","cog_sim_pred_hcap")
# out_dir <- "Analysis/Simulation Results/"
#
# r1 <- mergeSimResults(file_name = "ds_iteration_group_1.rds", out_dir = out_dir,
# dat_cols = dat_cols, res_cols = res_cols, iter_group = 50, niter = 20)
# saveRDS(r1,file=paste0(out_dir,"model_predicted_trajectory_simulations.rds"))
#
# rm(r1)
# prepare longitudinal data file
out_dir <- "~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/Simulation Results/10_15_24/"
dat_cols <- c("id","time","agebl_75","slope_true_qrtl","int_true_qrtl")
res_cols <- c("cog_true_pred_all","cog_sim_pred_all","cog_sim_pred_mmse",
"cog_sim_pred_tics","cog_sim_pred_hcap")
r1 <- readRDS(paste0(out_dir,"model_predicted_trajectory_simulations.rds"))
r1$cog_true_pred_mn <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_true_pred_se <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
# summary(r1$cog_true_mn)
r1$cog_sim_pred_all_mn <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_all_se <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_mmse_mn <- apply(r1[,names(r1)[grepl("sim_pred_mmse", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_mmse_se <- apply(r1[,names(r1)[grepl("sim_pred_mmse", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_tics_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics_se <- apply(r1[,names(r1)[grepl("sim_pred_tics", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_hcap_mn <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_hcap_se <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r2 <- r1 %>% dplyr::select(any_of(dat_cols),contains("_mn"),contains("_se")) %>%
mutate(age = agebl_75 + 75)
# summary(r2$age)
r3a <- r2 %>% dplyr::select(-c(cog_true_pred_se:cog_sim_pred_hcap_se)) %>%
pivot_longer(cols = cog_true_pred_mn:cog_sim_pred_hcap_mn,
names_to = "label", values_to = "cog") %>% mutate(
label = sub("_mn$","",label)
)
r3b <- r2 %>% dplyr::select(-c(cog_true_pred_mn:cog_sim_pred_hcap_mn)) %>%
pivot_longer(cols = cog_true_pred_se:cog_sim_pred_hcap_se,
names_to = "label", values_to = "cog_se") %>% mutate(
label = sub("_se$","",label)
)
r3 <- r3a %>% left_join((r3b %>%
dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>%
mutate(
cognitive_measure = factor(case_when(
label == "cog_true_pred" ~ "True Cognition",
label == "cog_sim_pred_all" ~ "MMSE+TICS+HCAP",
label == "cog_sim_pred_mmse" ~ "MMSE",
label == "cog_sim_pred_tics" ~ "TICS",
label == "cog_sim_pred_hcap" ~ "HCAP",
), levels = c("True Cognition","MMSE+TICS+HCAP","HCAP","TICS","MMSE"))
)
r4 <- r3 %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
mean_cog = mean(cog),
se_cog = mean(cog_se) / sqrt(1000)
) %>% mutate(
lower = mean_cog - (1.96 * se_cog),
upper = mean_cog + (1.96 * se_cog),
)
r2ta <- r2 %>% dplyr::select(id,time,agebl_75,slope_true_qrtl,
cog_sim_pred_tics_mn,cog_sim_pred_all_mn,cog_sim_pred_tics_se,cog_sim_pred_all_se) %>%
mutate(
cog_blend_mn = case_when(
time <= 5 ~ cog_sim_pred_tics_mn,
time > 5 ~ cog_sim_pred_all_mn
),
cog_blend_se = case_when(
time <= 5 ~ cog_sim_pred_tics_se,
time > 5 ~ cog_sim_pred_all_se
)
) %>% relocate(cog_blend_mn, .after = cog_sim_pred_all_mn)
r3taa <- r2ta %>% dplyr::select(-c(cog_sim_pred_tics_se:cog_blend_se,cog_sim_pred_all_mn)) %>%
pivot_longer(cols = c(cog_sim_pred_tics_mn,cog_blend_mn),
names_to = "label", values_to = "cog") %>% mutate(
label = sub("_mn$","",label)
)
r3tab <- r2ta %>% dplyr::select(-c(cog_sim_pred_tics_mn:cog_blend_mn,cog_sim_pred_all_se)) %>%
pivot_longer(cols = c(cog_sim_pred_tics_se,cog_blend_se),
names_to = "label", values_to = "cog_se") %>% mutate(
label = sub("_se$","",label)
)
r3ta <- r3taa %>% left_join((r3tab %>%
dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>%
mutate(
cognitive_measure = factor(case_when(
label == "cog_sim_pred_tics" ~ "TICS",
label == "cog_blend" ~ "Blended TICS,MMSE+TICS+HCAP"
), levels = c("Blended TICS,MMSE+TICS+HCAP","TICS"))
)
r4ta <- r3ta %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
mean_cog = mean(cog),
se_cog = mean(cog_se) / sqrt(1000)
) %>% mutate(
lower = mean_cog - (1.96 * se_cog),
upper = mean_cog + (1.96 * se_cog),
)
rm(r2,r3,r3a,r3b,r2ta,r3ta,r3taa,r3tab)
# cor(r1[,6:31])
### Repeat for iteration 1:100
r1 <- r1[,1:505]
r1$cog_true_pred_mn <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_true_pred_se <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
# summary(r1$cog_true_mn)
r1$cog_sim_pred_all_mn <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_all_se <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_mmse_mn <- apply(r1[,names(r1)[grepl("sim_pred_mmse", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_mmse_se <- apply(r1[,names(r1)[grepl("sim_pred_mmse", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_tics_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics_se <- apply(r1[,names(r1)[grepl("sim_pred_tics", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_hcap_mn <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_hcap_se <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r2 <- r1 %>% dplyr::select(any_of(dat_cols),contains("_mn"),contains("_se")) %>%
mutate(age = agebl_75 + 75)
# summary(r2$age)
r3a <- r2 %>% dplyr::select(-c(cog_true_pred_se:cog_sim_pred_hcap_se)) %>%
pivot_longer(cols = cog_true_pred_mn:cog_sim_pred_hcap_mn,
names_to = "label", values_to = "cog") %>% mutate(
label = sub("_mn$","",label)
)
r3b <- r2 %>% dplyr::select(-c(cog_true_pred_mn:cog_sim_pred_hcap_mn)) %>%
pivot_longer(cols = cog_true_pred_se:cog_sim_pred_hcap_se,
names_to = "label", values_to = "cog_se") %>% mutate(
label = sub("_se$","",label)
)
r3 <- r3a %>% left_join((r3b %>%
dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>%
mutate(
cognitive_measure = factor(case_when(
label == "cog_true_pred" ~ "True Cognition",
label == "cog_sim_pred_all" ~ "MMSE+TICS+HCAP",
label == "cog_sim_pred_mmse" ~ "MMSE",
label == "cog_sim_pred_tics" ~ "TICS",
label == "cog_sim_pred_hcap" ~ "HCAP",
), levels = c("True Cognition","MMSE+TICS+HCAP","HCAP","TICS","MMSE"))
)
r4a <- r3 %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
mean_cog = mean(cog),
se_cog = mean(cog_se) / sqrt(100)
) %>% mutate(
lower = mean_cog - (1.96 * se_cog),
upper = mean_cog + (1.96 * se_cog),
)
rm(r1,r2,r3,r3a,r3b)
### End data
pd <- position_dodge(width = 0.2) # move them .2 to the left and right
ggplot(data = r4, aes(x = time, y = mean_cog, color = slope_true_qrtl)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
facet_wrap(vars(cognitive_measure)) +
labs(caption = "True versus simulated cognitive trajectories (1000 Simulations) - All Cognition Measures")
ggplot(data = r4a, aes(x = time, y = mean_cog, color = slope_true_qrtl)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
facet_wrap(vars(cognitive_measure)) +
labs(caption = "True versus simulated cognitive trajectories (100 Simulations) - All Cognition Measures")
r4 <- r4 %>% mutate(n_sim = 1000)
r4a <- r4a %>% mutate(n_sim = 100)
r10 <- r4 %>% bind_rows(r4a) %>% mutate(
n_sim = factor(n_sim, levels = c(100,1000))
)
ggplot(data = r10, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = n_sim)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
facet_wrap(vars(cognitive_measure)) +
labs(caption = "True versus simulated cognitive trajectories by number of simulations - All Cognition Measures")
r5 <- r4 %>% filter(cognitive_measure %in% c("MMSE+TICS+HCAP","True Cognition"))
ggplot(data = r5, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = "True versus simulated cognitive change (1000 Simulations) - MMSE + TICS + HCAP")
r6 <- r4 %>% filter(cognitive_measure %in% c("HCAP","True Cognition"))
ggplot(data = r6, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = "True versus simulated cognitive change (1000 Simulations) - HCAP")
r7 <- r4 %>% filter(cognitive_measure %in% c("TICS","True Cognition"))
ggplot(data = r7, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = "True versus simulated cognitive change (1000 Simulations) - TICS")
r7a <- r4a %>% filter(cognitive_measure %in% c("TICS","True Cognition"))
ggplot(data = r7a, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = "True versus simulated cognitive change (100 Simulations) - TICS")
r8 <- r4 %>% filter(cognitive_measure %in% c("MMSE","True Cognition"))
ggplot(data = r8, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = "True versus simulated cognitive change (1000 Simulations) - MMSE")
r8a <- r4a %>% filter(cognitive_measure %in% c("MMSE","True Cognition"))
ggplot(data = r8a, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = "True versus simulated cognitive change (100 Simulations) - MMSE")
r9 <- r4 %>% filter(cognitive_measure %in% c("MMSE+TICS+HCAP","MMSE"))
ggplot(data = r9, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = "Simulated cognitive change (1000 Simulations) - MMSE + TICS + HCAP vs. MMSE")
ggplot(data = r4ta, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = str_wrap("TICS all assessments versus TICS (assessments bl-5) blended with MMSE+TICS+HCAP (assessments 6-10) (1000 Simulations)",width = 70))
rm(r4,r4a,r5,r6,r7,r7a,r8,r8a,r9,r10)
Simulated Cognitive Trajectories versus True Cognition Trajectories - TICS Variants
out_dir <- "~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/Simulation Results/tics/"
dat_cols <- c("id","time","agebl_75","slope_true_qrtl","int_true_qrtl")
res_cols <- c("cog_true_pred_tics3","cog_sim_pred_tics3","cog_sim_pred_tics10",
"cog_sim_pred_tics13","cog_sim_pred_tics16")
r1 <- readRDS(paste0(out_dir,"model_predicted_trajectory_simulations.rds"))
r1$cog_true_pred_mn <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_true_pred_se <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
# summary(r1$cog_true_mn)
r1$cog_sim_pred_tics3_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics3", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics3_se <- apply(r1[,names(r1)[grepl("sim_pred_tics3", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_tics10_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics10", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics10_se <- apply(r1[,names(r1)[grepl("sim_pred_tics10", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_tics13_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics13", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics13_se <- apply(r1[,names(r1)[grepl("sim_pred_tics13", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_tics16_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics16", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics16_se <- apply(r1[,names(r1)[grepl("sim_pred_tics16", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r2 <- r1 %>% dplyr::select(any_of(dat_cols),contains("_mn"),contains("_se")) %>%
mutate(age = agebl_75 + 75)
# summary(r2$age)
r3a <- r2 %>% dplyr::select(-c(cog_true_pred_se:cog_sim_pred_tics16_se)) %>%
pivot_longer(cols = cog_true_pred_mn:cog_sim_pred_tics16_mn,
names_to = "label", values_to = "cog") %>% mutate(
label = sub("_mn$","",label)
)
r3b <- r2 %>% dplyr::select(-c(cog_true_pred_mn:cog_sim_pred_tics16_mn)) %>%
pivot_longer(cols = cog_true_pred_se:cog_sim_pred_tics16_se,
names_to = "label", values_to = "cog_se") %>% mutate(
label = sub("_se$","",label)
)
r3 <- r3a %>% left_join((r3b %>%
dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>%
mutate(
cognitive_measure = factor(case_when(
label == "cog_true_pred" ~ "True Cognition",
label == "cog_sim_pred_tics3" ~ "TICS (wave 3)",
label == "cog_sim_pred_tics10" ~ "TICS (wave 3 + Number Series)",
label == "cog_sim_pred_tics13" ~ "TICS (wave 13)",
label == "cog_sim_pred_tics16" ~ "TICS (wave 16)"
), levels = c("True Cognition","TICS (wave 16)","TICS (wave 13)",
"TICS (wave 3 + Number Series)","TICS (wave 3)"))
)
r4 <- r3 %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
mean_cog = mean(cog),
se_cog = mean(cog_se) / sqrt(1000)
) %>% mutate(
lower = mean_cog - (1.96 * se_cog),
upper = mean_cog + (1.96 * se_cog),
)
r2ta <- r2 %>% dplyr::select(id,time,agebl_75,slope_true_qrtl,
cog_sim_pred_tics3_mn,cog_sim_pred_tics16_mn,cog_sim_pred_tics3_se,cog_sim_pred_tics16_se) %>%
mutate(
cog_blend_mn = case_when(
time <= 5 ~ cog_sim_pred_tics3_mn,
time > 5 ~ cog_sim_pred_tics16_mn
),
cog_blend_se = case_when(
time <= 5 ~ cog_sim_pred_tics3_se,
time > 5 ~ cog_sim_pred_tics16_se
)
) %>% relocate(cog_blend_mn, .after = cog_sim_pred_tics16_mn)
r3taa <- r2ta %>% dplyr::select(-c(cog_sim_pred_tics3_se:cog_blend_se,cog_sim_pred_tics16_mn)) %>%
pivot_longer(cols = c(cog_sim_pred_tics3_mn,cog_blend_mn),
names_to = "label", values_to = "cog") %>% mutate(
label = sub("_mn$","",label)
)
r3tab <- r2ta %>% dplyr::select(-c(cog_sim_pred_tics3_mn:cog_blend_mn,cog_sim_pred_tics16_se)) %>%
pivot_longer(cols = c(cog_sim_pred_tics3_se,cog_blend_se),
names_to = "label", values_to = "cog_se") %>% mutate(
label = sub("_se$","",label)
)
r3ta <- r3taa %>% left_join((r3tab %>%
dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>%
mutate(
cognitive_measure = factor(case_when(
label == "cog_sim_pred_tics3" ~ "TICS (wave 3)",
label == "cog_blend" ~ "Blended TICS (wave 3),TICS (wave 16)"
), levels = c("Blended TICS (wave 3),TICS (wave 16)","TICS (wave 3)"))
)
r4ta <- r3ta %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
mean_cog = mean(cog),
se_cog = mean(cog_se) / sqrt(1000)
) %>% mutate(
lower = mean_cog - (1.96 * se_cog),
upper = mean_cog + (1.96 * se_cog),
)
rm(r2,r3,r3a,r3b,r2ta,r3ta,r3taa,r3tab)
ggplot(data = r4, aes(x = time, y = mean_cog, color = slope_true_qrtl)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
facet_wrap(vars(cognitive_measure)) +
labs(caption = "True versus simulated cognitive trajectories (1000 Simulations) - All TICS Measures")
ggplot(data = r4ta, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = str_wrap("TICS (wave 3) all assessments versus TICS (wave 3 - assessments bl-5) blended with TICS (wave 16 - assessments 6-10) (1000 Simulations)",width = 70))
